UW Summer Institutes - Introduction to Metagenomic Data Analysis - Lab 2

Getting ready

library(dada2); packageVersion("dada2")

[1] ‘1.0.3’

library(phyloseq); packageVersion("phyloseq")

[1] ‘1.16.2’

library(stringr); packageVersion("stringr")

[1] ‘1.0.0’

library(ggplot2); packageVersion("ggplot2")

[1] ‘2.1.0’

library(dplyr); packageVersion("dplyr")

[1] ‘0.5.0’

Load data from Nature 488, pp. 621-626: STAT.RData

path <- "../data/STAT.rdata"
load(path)
phy

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 6547 taxa and 96 samples ]
sample_data() Sample Data:       [ 96 samples by 6 sample variables ]
tax_table()   Taxonomy Table:    [ 6547 taxa by 9 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 6547 tips and 6321 internal nodes ]

head(phenotypes)

     GIP Insulin Leptin IGF1    BMD mFat pFat
C1  8.61       1   2320 1250 0.0492  4.0 19.7
C2 28.60     137   1040 2180 0.0487  3.7 18.6
C3 19.20     543   1670 2160 0.0510  3.6 19.5
C4 38.30     606    972 2070 0.0465  4.0 21.3
C5 13.20       1   2210 1240 0.0547  4.1 21.6
C6 36.80     696   2330 1830 0.0512  3.7 17.7

tail(phenotypes)

      GIP Insulin Leptin IGF1    BMD mFat pFat
VP5  49.5     418   2130 1810 0.0484  4.2 21.9
VP6  49.9     271   1740 1610 0.0476  4.7 23.9
VP7  48.6     224   1710 1520 0.0482  5.5 26.1
VP8  37.5       1    795 2020 0.0504  4.0 19.9
VP9  41.7    1200   4310 1390 0.0484  4.8 27.6
VP10 49.1       1    917 1970 0.0509  5.3 25.5

head(sample_data(phy))

Sample Data:        [6 samples by 6 sample variables]:
          X.SampleID BarcodeSequence LinkerPrimerSequence Treatment Location
cecal_C1    cecal_C1              NA                   NA         C    cecal
cecal_C10  cecal_C10              NA                   NA         C    cecal
cecal_C2    cecal_C2              NA                   NA         C    cecal
cecal_C3    cecal_C3              NA                   NA         C    cecal
cecal_C4    cecal_C4              NA                   NA         C    cecal
cecal_C5    cecal_C5              NA                   NA         C    cecal
                  Description
cecal_C1  missing_description
cecal_C10 missing_description
cecal_C2  missing_description
cecal_C3  missing_description
cecal_C4  missing_description
cecal_C5  missing_description

levels(sample_data(phy)$Treatment)

[1] "C"  "P"  "T"  "V"  "VP"

Merge phenotypes variable with the phyloseq data

unique(sample_names(phy))

 [1] "cecal_C1"   "cecal_C10"  "cecal_C2"   "cecal_C3"   "cecal_C4"   "cecal_C5"  
 [7] "cecal_C6"   "cecal_C7"   "cecal_C8"   "cecal_C9"   "cecal_P1"   "cecal_P10" 
[13] "cecal_P2"   "cecal_P3"   "cecal_P4"   "cecal_P5"   "cecal_P6"   "cecal_P7"  
[19] "cecal_P8"   "cecal_P9"   "cecal_T1"   "cecal_T10"  "cecal_T2"   "cecal_T3"  
[25] "cecal_T4"   "cecal_T5"   "cecal_T6"   "cecal_T7"   "cecal_T8"   "cecal_T9"  
[31] "cecal_V1"   "cecal_V10"  "cecal_V2"   "cecal_V3"   "cecal_V4"   "cecal_V5"  
[37] "cecal_V6"   "cecal_V7"   "cecal_V8"   "cecal_V9"   "cecal_VP1"  "cecal_VP10"
[43] "cecal_VP2"  "cecal_VP3"  "cecal_VP4"  "cecal_VP5"  "cecal_VP6"  "cecal_VP7" 
[49] "cecal_VP8"  "cecal_VP9"  "fecal_C1"   "fecal_C10"  "fecal_C2"   "fecal_C3"  
[55] "fecal_C4"   "fecal_C5"   "fecal_C6"   "fecal_C7"   "fecal_C8"   "fecal_C9"  
[61] "fecal_P1"   "fecal_P10"  "fecal_P2"   "fecal_P3"   "fecal_P4"   "fecal_P6"  
[67] "fecal_P7"   "fecal_P8"   "fecal_P9"   "fecal_T1"   "fecal_T10"  "fecal_T2"  
[73] "fecal_T3"   "fecal_T4"   "fecal_T6"   "fecal_T7"   "fecal_T8"   "fecal_T9"  
[79] "fecal_V1"   "fecal_V10"  "fecal_V2"   "fecal_V3"   "fecal_V4"   "fecal_V5"  
[85] "fecal_V6"   "fecal_V7"   "fecal_V8"   "fecal_V9"   "fecal_VP1"  "fecal_VP10"
[91] "fecal_VP2"  "fecal_VP3"  "fecal_VP4"  "fecal_VP7"  "fecal_VP8"  "fecal_VP9" 

pheno1 <- phenotypes
pheno2 <- phenotypes
rownames(pheno1) <- paste('cecal', rownames(pheno1), sep = "_")
rownames(pheno2) <- paste('fecal', rownames(pheno2), sep = "_")
pheno <- rbind(pheno1, pheno2)
head(pheno)

           GIP Insulin Leptin IGF1    BMD mFat pFat
cecal_C1  8.61       1   2320 1250 0.0492  4.0 19.7
cecal_C2 28.60     137   1040 2180 0.0487  3.7 18.6
cecal_C3 19.20     543   1670 2160 0.0510  3.6 19.5
cecal_C4 38.30     606    972 2070 0.0465  4.0 21.3
cecal_C5 13.20       1   2210 1240 0.0547  4.1 21.6
cecal_C6 36.80     696   2330 1830 0.0512  3.7 17.7

tail(pheno)

            GIP Insulin Leptin IGF1    BMD mFat pFat
fecal_VP5  49.5     418   2130 1810 0.0484  4.2 21.9
fecal_VP6  49.9     271   1740 1610 0.0476  4.7 23.9
fecal_VP7  48.6     224   1710 1520 0.0482  5.5 26.1
fecal_VP8  37.5       1    795 2020 0.0504  4.0 19.9
fecal_VP9  41.7    1200   4310 1390 0.0484  4.8 27.6
fecal_VP10 49.1       1    917 1970 0.0509  5.3 25.5

# phenotypes_fixed <- phenotypes
# rownames(phenotypes_fixed) <- str_replace_all(rownames(phenotypes), 
#                                               "C", "cecal_C")
# rownames(phenotypes_fixed) <- str_replace_all(rownames(phenotypes_fixed), 
#                                               "C", "cecal_C")

phy2 <- merge_phyloseq(phy, sample_data(pheno))
phy

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 6547 taxa and 96 samples ]
sample_data() Sample Data:       [ 96 samples by 6 sample variables ]
tax_table()   Taxonomy Table:    [ 6547 taxa by 9 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 6547 tips and 6321 internal nodes ]

phy2

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 6547 taxa and 96 samples ]
sample_data() Sample Data:       [ 96 samples by 13 sample variables ]
tax_table()   Taxonomy Table:    [ 6547 taxa by 9 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 6547 tips and 6321 internal nodes ]

head(sample_data(phy2))

Sample Data:        [6 samples by 13 sample variables]:
          X.SampleID BarcodeSequence LinkerPrimerSequence Treatment Location
cecal_C1    cecal_C1              NA                   NA         C    cecal
cecal_C10  cecal_C10              NA                   NA         C    cecal
cecal_C2    cecal_C2              NA                   NA         C    cecal
cecal_C3    cecal_C3              NA                   NA         C    cecal
cecal_C4    cecal_C4              NA                   NA         C    cecal
cecal_C5    cecal_C5              NA                   NA         C    cecal
                  Description   GIP Insulin Leptin IGF1    BMD mFat pFat
cecal_C1  missing_description  8.61       1   2320 1250 0.0492  4.0 19.7
cecal_C10 missing_description 27.70     516   2680 3720 0.0479  4.0 19.8
cecal_C2  missing_description 28.60     137   1040 2180 0.0487  3.7 18.6
cecal_C3  missing_description 19.20     543   1670 2160 0.0510  3.6 19.5
cecal_C4  missing_description 38.30     606    972 2070 0.0465  4.0 21.3
cecal_C5  missing_description 13.20       1   2210 1240 0.0547  4.1 21.6

Estimate observed species with Chao1 and Shannon diversity using estimate_richness() function

alpha_meas <- c("Observed", "Chao1", "Shannon")
alpha_div <- estimate_richness(phy2, measures = alpha_meas)
glimpse(alpha_div)

Observations: 96
Variables: 4
$ Observed <dbl> 662, 592, 617, 731, 601, 698, 579, 529, 462, 617, 484, 593, ...
$ Chao1    <dbl> 1650.4286, 1508.4531, 1216.3939, 1746.2188, 1452.0685, 1434....
$ se.chao1 <dbl> 148.60209, 153.81056, 91.48394, 144.87205, 137.36109, 104.96...
$ Shannon  <dbl> 4.836219, 4.814829, 4.775428, 4.918447, 4.636482, 4.944320, ...

Subsetting

phy3 <- subset_samples(phy2, sample_sums(phy2) > 3500)
phy3

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 6547 taxa and 93 samples ]
sample_data() Sample Data:       [ 93 samples by 13 sample variables ]
tax_table()   Taxonomy Table:    [ 6547 taxa by 9 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 6547 tips and 6321 internal nodes ]

alpha_div_sub <- estimate_richness(phy3, measures = alpha_meas)

Estimate rarified diversities in #3 averaged over 20 replicates and compare the two estimates

set.seed(0)
alpha_div_rare <- matrix(0, ncol = ncol(alpha_div), nrow = nrow(alpha_div))
for (i in 1:20) {
    alpha_div_rare <- alpha_div_rare +
        estimate_richness(
            rarefy_even_depth(phy3, sample.size = 3500, 
                              verbose = F, trimOTUs = F),
            measures = alpha_meas
        )
}
alpha_div_rare <- alpha_div_rare / 20
glimpse(alpha_div_rare)

Observations: 93
Variables: 4
$ Observed <dbl> 467.70, 409.80, 443.30, 501.20, 398.80, 481.90, 376.10, 354....
$ Chao1    <dbl> 643.5768, 550.6086, 598.5164, 705.6298, 563.2645, 668.5343, ...
$ se.chao1 <dbl> 33.85371, 29.96853, 31.54622, 37.52648, 34.22859, 34.91094, ...
$ Shannon  <dbl> 4.732475, 4.723733, 4.675471, 4.806043, 4.540449, 4.826177, ...

plot(alpha_div$Observed, alpha_div_rare$Observed)

Error in xy.coords(x, y, xlabel, ylabel, log) : 
  'x' and 'y' lengths differ

Summarize the data at the order level

order_phy <- tax_glom(phy3, taxrank = "Order")
glimpse(order_phy)

Formal class 'phyloseq' [package "phyloseq"] with 5 slots
  ..@ otu_table:Formal class 'otu_table' [package "phyloseq"] with 2 slots
  ..@ tax_table:Formal class 'taxonomyTable' [package "phyloseq"] with 1 slot
  ..@ sam_data :'data.frame':   93 obs. of  13 variables:
Formal class 'sample_data' [package "phyloseq"] with 4 slots
  ..@ phy_tree :List of 5
  .. ..$ edge       : int [1:39, 1:2] 22 22 23 24 24 25 26 26 27 27 ...
  .. ..$ Nnode      : int 19
  .. ..$ tip.label  : chr [1:21] "1948" "1931" "4911" "4220" ...
  .. ..$ edge.length: num [1:39] 0.2476 0.2294 0.1725 0.4988 0.0353 ...
  .. ..$ node.label : chr [1:19] "" "0.783" "0.768" "0.334" ...
  .. ..- attr(*, "class")= chr "phylo"
  .. ..- attr(*, "order")= chr "cladewise"
  ..@ refseq   : NULL

Normalize the data using CLR, relative abundance, and DESeq2

order_rel <- transform_sample_counts(order_phy,
                                     function(x) x / sum(x))

order_clr <- transform_sample_counts(order_phy,
                                     function(x) {
                                         y = log(x + 1)
                                         y / sum(y)
                                     })

glimpse(order_deseq)

Formal class 'phyloseq' [package "phyloseq"] with 5 slots
  ..@ otu_table:Formal class 'otu_table' [package "phyloseq"] with 2 slots
  ..@ tax_table:Formal class 'taxonomyTable' [package "phyloseq"] with 1 slot
  ..@ sam_data :'data.frame':   93 obs. of  13 variables:
Formal class 'sample_data' [package "phyloseq"] with 4 slots
  ..@ phy_tree :List of 5
  .. ..$ edge       : int [1:39, 1:2] 22 22 23 24 24 25 26 26 27 27 ...
  .. ..$ Nnode      : int 19
  .. ..$ tip.label  : chr [1:21] "1948" "1931" "4911" "4220" ...
  .. ..$ edge.length: num [1:39] 0.2476 0.2294 0.1725 0.4988 0.0353 ...
  .. ..$ node.label : chr [1:19] "" "0.783" "0.768" "0.334" ...
  .. ..- attr(*, "class")= chr "phylo"
  .. ..- attr(*, "order")= chr "cladewise"
  ..@ refseq   : NULL

Compute appropriate univariate tests on normalized data

Location = sample_data(order_rel)$Location
rel.p = apply(otu_table(order_rel),1, function(x) wilcox.test(c(x)~Location)$p.value)
clr.p = apply(otu_table(order_clr),1, function(x) t.test(c(x)~Location)$p.value)
deseq.p = apply(otu_table(order_deseq),1, function(x) t.test(log10(c(x))~Location)$p.value)

univ.location.res = data.frame(taxa = tax_table(order_clr)[,"Order"], rel.p, clr.p, deseq.p)
head(univ.location.res)

Adjust the p-values using False Discovery Rate

univ.location.res$rel.fdr = p.adjust(univ.location.res$rel.p, method="fdr")
univ.location.res$clr.fdr = p.adjust(univ.location.res$clr.p, method="fdr")
univ.location.res$deseq.fdr = p.adjust(univ.location.res$deseq.p, method="fdr")
colSums(univ.location.res<0.05)

Warning in Ops.factor(left, right) : ‘<’ not meaningful for factors

    Order     rel.p     clr.p   deseq.p   rel.fdr   clr.fdr deseq.fdr 
       NA        10         8         9         9         7         7 

univ.location.res

                    Order        rel.p        clr.p      deseq.p      rel.fdr
1948        Bacteroidales 1.453797e-11 9.983697e-01 4.902762e-01 8.177357e-11
1931 "Erysipelotrichales" 8.437764e-01 2.076959e-01 5.187488e-01 8.859653e-01
4911    Anaeroplasmatales 9.721986e-10 2.249341e-15 5.735335e-14 4.083234e-09
4220        Clostridiales 1.186926e-19 1.894648e-32 3.651509e-18 2.492544e-18
1026           Bacillales 8.032710e-08 1.301173e-06 4.190743e-06 2.811448e-07
1673    "Lactobacillales" 6.834879e-16 4.321846e-19 1.071339e-16 7.176623e-15
4789     Actinobacteridae 1.132115e-01 7.601720e-02 4.509191e-02 1.828801e-01
5995      Burkholderiales 7.245135e-02 8.371138e-02 1.021904e-01 1.269279e-01
5975         Neisseriales 3.437768e-01 3.224313e-01 1.034844e-01 4.010729e-01
5242     Sphingomonadales 7.253024e-02 9.877532e-02 2.072876e-01 1.269279e-01
1453      Rhodobacterales 3.437768e-01 3.224313e-01 1.066468e-01 4.010729e-01
455           Rhizobiales 3.118829e-01 3.227785e-01 8.148514e-01 4.010729e-01
6648         Myxococcales 3.118829e-01 3.227785e-01 8.165448e-01 4.010729e-01
3405    Enterobacteriales 2.819788e-02 3.464769e-02 3.530636e-02 5.921555e-02
2012       Pasteurellales 9.633972e-01 9.332040e-01 6.386468e-01 9.633972e-01
1847      Pseudomonadales 1.239806e-03 4.114033e-03 4.024734e-03 3.254492e-03
2172          Chloroplast 8.788168e-05 3.013895e-04 2.380779e-04 2.636450e-04
3806     Flavobacteriales 1.461646e-01 1.886147e-01 3.350350e-01 2.192469e-01
4609      Mycoplasmatales 5.427940e-01 2.501740e-01 2.868498e-01 5.999302e-01
3348   Desulfovibrionales 1.954624e-02 9.799033e-01 2.337175e-01 4.560790e-02
2626      Coriobacteridae 1.557592e-11 1.459745e-12 3.481561e-12 8.177357e-11
          clr.fdr    deseq.fdr
1948 9.983697e-01 6.052069e-01
1931 3.355087e-01 6.052069e-01
4911 1.574539e-14 4.014735e-13
4220 3.978760e-31 7.668168e-17
1026 5.464927e-06 1.760112e-05
1673 4.537938e-18 1.124905e-15
4789 1.757939e-01 1.052145e-01
5995 1.757939e-01 1.866319e-01
5975 3.765749e-01 1.866319e-01
5242 1.885711e-01 3.348492e-01
1453 3.765749e-01 1.866319e-01
455  3.765749e-01 8.165448e-01
6648 3.765749e-01 8.165448e-01
3405 9.095018e-02 9.267919e-02
2012 9.983697e-01 7.058727e-01
1847 1.234210e-02 1.207420e-02
2172 1.054863e-03 8.332725e-04
3806 3.300758e-01 4.397334e-01
4609 3.752611e-01 4.015897e-01
3348 9.983697e-01 3.505762e-01
2626 7.663660e-12 1.827820e-11

Plot the abundances using plot_heatmap(), plot_bar(), or produce box and whisker plots

plot_heatmap(order_rel, taxa.label = "Order", sample.order = "Location")+
  facet_grid(.~Location, scales = "free_x")

plot_heatmap(order_clr, taxa.label = "Order", sample.order = "Location")+
  facet_grid(.~Location, scales = "free_x")

plot_heatmap(order_deseq, taxa.label = "Order", sample.order = "Location")+
  facet_grid(.~Location, scales = "free_x")

plot_bar(order_rel, fill="Location") + 
  facet_wrap(facets = "Order", ncol=5, nrow=5, scales = "free_y")

plot_bar(order_clr, fill="Location") + 
  facet_wrap(facets = "Order", ncol=5, nrow=5, scales = "free_y")

plot_bar(order_deseq, fill="Location") + 
  facet_wrap(facets = "Order", ncol=5, nrow=5, scales = "free_y")

ggplot(psmelt(order_rel), aes(x=Location, y=Abundance)) + 
  geom_boxplot(aes(fill=Location)) + geom_jitter() + 
  facet_wrap(facets="Order", ncol=5, nrow=5, scales="free_y")

ggplot(psmelt(order_clr), aes(x=Location, y=Abundance)) + 
  geom_boxplot(aes(fill=Location)) + geom_jitter() + 
  facet_wrap(facets="Order", ncol=5, nrow=5, scales="free_y")

ggplot(psmelt(order_deseq), aes(x=Location, y=Abundance)) + 
  geom_boxplot(aes(fill=Location)) + geom_jitter() + 
  facet_wrap(facets="Order", ncol=5, nrow=5, scales="free_y")

Subset the order level data to only fecal samples

fecal.rel = subset_samples(order_rel, Location = "fecal")
fecal.rel.subs = subset_taxa(fecal.rel, rowMeans(otu_table(fecal.rel))> 0.001)

Filter out taxa with low abundance (mean abundance < 0.1%)

fecal.rel.subs = subset_taxa(fecal.rel, rowMeans(otu_table(fecal.rel))> 0.001)

Perform univariate analysis of the fecal subset with respect to the Treatment variable